Pair dynamics in a glass forming binary mixture: Simulations and theory 
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We have carried out molecular dynamics simulations to understand the dynamics of a tagged pair 
of atoms in a strongly non-ideal glass-forming binary Lennard- Jones mixture. Here atom B is smaller 
than atom A (cjbb = 0.88a aa, where oaa is the molecular diameter of the A particles) and the 
AB interaction is stronger than that given by Lorentz-Berthelot mixing rule {cab = IJatAA, where 
caa is the interaction energy strength between the A particles). The generalized time-dependent 
pair distribution function is calculated separately for the three pairs (AA, BB and AB). The three 
pairs are found to behave differently. The relative diffusion constants are found to vary in the order 
D BB > D AB > Dft A , with Dr B ~ 2Dr A , showing the importance of the hopping process (B 
hops much more than A). We introduce a non-Gaussian parameter (a P (t)) to monitor the relative 
motion of a pair of atoms, and evaluate it for all the three pairs, with initial separations chosen to be 
at the first peak of the corresponding partial radial distribution functions. At intermediate times, 
significant deviation from the Gaussian behavior of the pair distribution functions is observed, with 
different degree for the three pairs. A simple mean-field (MF) model, proposed originally by Haan 
[Phys. Rev. A 20, 2516 (1979)] for one component liquid, is applied to the case of binary mixture, 
and compared with the simulation results. While the MF model successfully describe the dynamics 
of the AA and AB pair, the agreement for the BB pair is less satisfactory. This is attributed to the 
large scale anharmonic motions of the B particles in a weak effective potential. Dynamics of next 
nearest neighbor pairs are also investigated. 



I. INTRODUCTION 



In dense fluids, there are many interaction-induced 
phenomena that can be interpreted in terms of the dy- 
namics of the pairs of atoms[jl ||, [|. For example, nu- 
clear overheusser effect studies the relative motion of the 
atoms. In addition, an understanding of pair dynam- 
ics can be of great importance in the studies of rate of 
various diffusion controlled chemical reactions in dense 
fluids d |. Both the theoretical analysis J& |, @, |, § [u], 
[Tl| and computer simulation^, |], [r], [l^] studies have 
been carried out extensively to study the dynamics of a 
pair of atoms in an one component liquid. Surprisingly, 
however, we are not aware of any explicit study of the 
dynamics of atomic pairs in binary mixtures whose dy- 
namics generally show strong nonmonotonic composition 
dependence @ Q. 

The study of the electronic spectroscopy of dilute chro- 
mophores ('solutes') in fluids ('solvents') is a useful tool 
for obtaining the information about the structure and 
dynamics of the solvents in the vicinity of the solute. In 
an attempt to provide a microscopic foundation of the 
Kubo's stochastic theory of the line shape, Skinner and 
coworkers |ll| have recently developed a molecular theory 
for the absorption and emission line shapes and ultrafast 
solvation dynamics of a dilute nonpolar solute in nonpo- 
lar fluids. Due to the motion of the solvent molecules rel- 
ative to the chromophore, the chromophore's transition 
frequency generally fluctuates in time. Thus the nature 
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of the spectral line shape provide a useful information 
about the details of the dynamics of the solvent relative 
to the solute. An approximate treatment of the solvent 
dynamics allowed the theory to express the transition 
frequency fluctuation time correlation functions (related 
to the expressions for the absorption and emission line 
shapes) solely in terms of the two-body solute-solvent 
time-dependent conditional pair distribution function. 
Many other applications of pair dynamics have been dis- 
cussed by a number of authors |6[ [To[ [Tl| [l6||. 

The dynamics of a liquid below its freezing temper- 
ature, that is, in a supercooled state, is far more com- 
plex than what one would expect from a naive extrap- 
olation of their high-temperature behavior. One of the 
most challenging problems in the dynamics of a super- 
cooled liquid is to understand quantitatively the origin 
of the non-exponential relaxation exhibited by various 
dynamical response functions and the extraordinary vis- 
cous slow-down within a narrow temperature range as 
one approaches the glass transition temperature from 
above |l7| ^8|. Many experimental studies jl9], as 
well as computer simulations 23, Q have been 
performed to shade light on the underlying microscopic 
mechanism involved in supercooled liquids. These stud- 
ies have revealed the evidence of the presence of dis- 
tinct relaxing domains (spatial heterogeneity) which is 
thought to be responsible for the non-exponential relax- 
ations in deeply supercooled liquids. Molecular motions 
in strongly supercooled liquid involves highly collective 
movement of several molecules^, [2^, |27], ^8|. Fur- 
thermore, the correlated jump motions become the dom- 
inant diffusive mode (2^, |2^]. The observed heterogeneity 
of the relaxations in a deeply supercooled liquid is found 
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to be connected to the collective hopping of groups of 
particles |30|. 

The occurrence of increasingly heterogeneous dynam- 
ics in supercooled liquids, however, has been investigated 
solely in terms of single particle dynamics. The study of 
the dynamics of pair of atoms which involve higher order 
(two-body) correlations thus can provide much broader 
insight into the anomalous dynamics of supercooled liq- 
uids. In this work, we have carried out molecular dy- 
namics simulations in a strongly non-ideal glass form- 
ing binary mixture (commonly known as Kob-Andersen 
model |pl| ) to study the relaxation mechanism in terms 
of the pair dynamics. The main purpose of the present 
study is to explore the dynamics in a more collective 
sense by following the relative motion of three different 
type (AA, BB, and AB) of nearest neighbor and next 
nearest neighbor pair of atoms. These three pairs are 
found to behave differently. The simulation results show 
a clear signature of hooping motion in all the three pairs. 
We have also performed simple mean-field (MF) model 
(as introduced by Haanj(| for one component liquid) cal- 
culations to obtain the time dependent conditional pair 
distribution functions. 

The organization of the rest of the paper is as follows. 
In Sec. II, we describe the details of the simulation and 
the model system used in this study. The simulation 
results are presented and discussed in Sec. III. In Sec. 
IV, we have presented a mean-field model calculations 
for pair dynamics in binary mixture and the comparison 
is made with the simulation results. Finally, we end with 
a few concluding remarks in Sec. V. 



order to lower the computational burden the potential 
has been truncated with a cutoff radius of 2.5a aa- All 
the quantities in this study are given in reduced units, 
that is, length in units of <jaa, temperature T in units 
of eAA/ks, pressure P in units of caa/o~aa- The corre- 
sponding microscopic time scale is r = \J ma\ A j ^aa- 

All simulations in the NPT ensemble were performed 
using the Nose-Hoover- Andersen method^] where the 
external reduced temperature is held fixed at T* = 1.0. 
The external reduced pressure has been kept fixed at 
P* = 20.0. The reduced average density p*(= P^aa) 01 
the system corresponding to this thermodynamic state 
point is 1.32. Throughout the course of the simulations, 
the barostat and the system's degrees of freedom are 
coupled to an independent Nose-Hoover chain f^if (NHC) 
of thermostats, each of length 5. The extended system 
equations of motion are integrated using the reversible 
integrator method|Q with a time step of 0.002. The 
higher order multiple time step method has been em- 
ployed in the NHC evolution operator which lead to sta- 
ble energy conservation for non-Hamiltonian dynamical 
systems [f35|. The extended system time scale parameter 
used in the calculations was taken to be 1.15 for both the 
barostat and thermostats. 

The systems were equilibrated for 2 x 10 6 time steps 
and simulations were carried out for another 10 7 pro- 
duction steps following equilibration, during which the 
quantities of interest are calculated. 



III. SIMULATION RESULTS AND DISCUSSION 



II. SYSTEM AND SIMULATION DETAILS 

We performed a series of equilibrium isothermal- 
isobaric ensemble (N P T) molecular dynamics (MD) sim- 
ulation of a strongly non-ideal well-known glass forming 
binary mixture in three dimensions. The binary system 
studied here contains a total of N = 1000 particles con- 
sisting of two species of particles, A and B with Na = 800 
and Nb — 200 number of A and B particles, respectively. 
Thus, the mixture consists of 80% of A particles and 20% 
of B particles. The interaction between any two particles 
is modeled by means of shifted force Lennard- Jones (LJ) 
pair potential EJ, where the standard LJ is given by 
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where i and j denote two different particles (A and B). 
The potential parameters are as follows: eaa = 1.0, 
cjaa = 1-0, tBB = 0.5, obb = 0.88, e A B = 1-5, 
and <7i2 = 0.8. The mass of the two species are same 
(in a — ttib = m). Note that in this model system the AB 
interaction (cab) is much stronger than both of their re- 
spective pure counterparts and a ab is smaller than what 
is expected from the Lorentz-Berthelot mixing rules. In 



The three partial radial distribution functions, qaa(t), 
9ab(t) and (7ss(r) obtained from simulations are plot- 
ted in figure 1. Due to the strong mutual interaction, 
the AB correlation is obviously the strongest among the 
three pairs. The splitting of the second peak of both 
9aa(t) and qab{t) is a characteristic signature of dense 
random packing. The structure of 3ss(r) is interest- 
ingly different. It has an insignificant first peak which 
originates from the weak interaction between the B type 
particles. The second peak of Qbb (t) is higher than that 
of the first peak signifying that the predominant BB cor- 
relation takes place at the second coordination shell. The 
occurrence of the splitted second peak is clearly observed 
here also. 

In the present study, the central quantity of interest is 
the time- dependent pair distribution function (TDPDF) 
(first introduced by Oppenheim and Bloomjp in the the- 
ory of nuclear spin relaxation in fluids), gi (r , r; t) which 
is the conditional probability that two particles are sep- 
arated by r at time t if that pair were separated by r Q 
at time t = 0. Thus, the TDPDF measures the relative 
motion of a pair of atoms. For an isotropic fluid, the 
TDPDF depends only on the magnitudes of r, r D and 

where 9 is the angle between r and r Q . In computer 
simulation, one can readily evaluate separately the radial 
and orientational features of the relative motion. In the 
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following two subsections, we present, respectively, the 
results obtained for the time evolution of the radial part 
g2,rad(r ,r;t) and the angular part g2,ang(r ,9;t) of the 
TDPDF for the three different pairs (AA, BB and AB). 



A. Radial part of the TDPDF, g 2 ,rad(ro, r; t) 

In figure 2 we plot the g2,rad(r 0l r;t) for the AA 
pair with the initial separation r a corresponds to the 
first maximum of the partial radial distribution func- 
tion gAA {r) (that is, pairs which are nearest neighbor) 
at four different times. While at short time (figure 
2(a)) the distribution function has a single peak struc- 
ture as expected, with increase in time it reaches slowly 
to its asymptotic limit where additional peaks develop 
at larger relative separations (see figures 2(b)-2(d)). The 
microscopic details of the underlying diffusive process (by 
which it approach to the asymptotic structure) can be ob- 
tained by following the trajectory of the relative motions. 
Figure 3 displays the projections onto an x-y plane of the 
trajectory of a typical AA pair for the nearest neighbor 
A atoms over a time interval At = 500r. The motion of 
the AA pair is shown to be relatively localized for many 
time steps and then move significant distances only dur- 
ing quick, rare cage rearrangements. This is a clear evi- 
dence that the jump motions are the dominant diffusive 
mode by which the separation between pairs of atoms 
evolve in time. 

The behavior of the distribution function g2,rad(r , r; t) 
for the AB pairs (where the interaction being the 
strongest) is plotted in figure 4 at four different times. 
The distribution function shows the same qualitative be- 
havior as we observed in the case of AA correlation (fig- 
ure 2). When compared to the AA correlation function 
within the same time scale, the decay of the correlation 
function is found to be faster despite the much strong AB 
interaction. This must be attributed to the difference in 
size of the two types of particles. As the B particles are 
smaller in size than the A particles, they are more mo- 
bile. In addition, the AB interaction is such that AB 
repulsion is felt at relatively small distances (<tab = 0.8 
instead of 0.94 according to the Lorentz-Berthelot rules). 
Consequently, the B particles are more prone to make 
jumps than the A particles (as observed earlier by Kob 
and Andersen|pTf). 

The nature of the relative motion of a typical AB pair 
is illustrated in figure 5(a), which display the trajectory 
of a typical AB pair (in the x-y plane) that were ini- 
tially at the nearest neighbor (first peak of gAB(r)). The 
elapsed time is At = 500r. The dynamics of the rela- 
tive motion is again dominated by hooping, the AB pair 
remain trapped at their initial separation over hundreds 
of time steps, before jumping to neighboring sites where 
they again become localized. Further, the jump motion 
is more frequent for the AB pair than that for the AA 
pair. The individual trajectory of the A and B parti- 
cles of the same AB pair within the same time window is 



shown in figures 5(b) and 5(c), respectively. While both 
the A and B particles hop, B particles move faster and 
the effect of caging is weaker (than the A particles) due 
to its smaller size. In this time window, the net displace- 
ment of the AB pair in the x-y plane is found to be quite 
large as shown in figure 6 and mainly determined by the 
displacement of the B particle. 

In figure 7 we show g2,rad(r , r; t) for the BB pair ini- 
tially separated at the first peak of gs b i r ) at four differ- 
ent times. Due to weak interaction among B particles, 
one expects that the B atoms in the BB pair will fast 
lose the memory of their initial separation. This is in- 
deed the case for the BB pair shown in figure 7. Once 
again the jump dynamics is clearly seen in the trajectory 
of a typical BB pair projected in the x-y plane (figure 8). 

We now consider the case where the separation of the 
initial pairs corresponds to the second peak of their re- 
spective partial radial distribution functions of figure 1 
(that is, pairs which are next nearest neighbor). The 
distribution function for the AA pair is plotted in fig- 
ure 9. It shows a qualitative different behavior because 
the peak at the nearest neighbor separation develops in 
a relatively short time. Here also the motion of the pairs 
are found to be mostly discontinuous in nature, thus mo- 
tion from second to first nearest neighbor occurs mostly 
by hopping. In figure 10 we plot the similar distribution 
function of the AB pair. Since the AB interaction is the 
strongest, the height of the first peak grows faster than 
that for the AA pair (compare figures 9(b) and 10(b)). 
Next, in figure 11 we plot the distribution function for the 
BB pair. Contrary to the AA and AB pairs, BB pairs 
tend to retain their initial separation for a relatively long 
time compared to the nearest neighbor pair. This can be 
understood from the predominant BB correlations at the 
second coordination shell. 



B. Angular part of the TDPDF, g2,ang{r o ,0;t) 

In this subsection we present the angular distribution 
function g2,ang(j' ,9;t) for the three different pairs (AA, 
BB and AB). The initial separation r Q for the three pairs 
corresponds to the first peak of the respective partial 
radial distribution functions (figure 1). 

In figure 12(a) we show the angular distribution 
g2,ang(r , Q\ t) for the AA pair. We calculate the angular 
distribution with respect to the initial separation vector 
r Q and irrespective of the value of the separation at time 
t. The distribution which is a 5-function at t = 0, spreads 
more and more with time and eventually it reaches to 
a uniform distribution with zero slope. When we com- 
pare to the distribution corresponding to the AB pair 
as shown in figure 12(b), we find that the approach to 
the uniform value is faster in the case of AB pair. This 
can be understood again in terms of the mobility of the 
B particles which is more compared to the A particles. 
In figure 12(c) we show how the distribution for the BB 
pair changes with time. The relaxation is seen to be rel- 
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atively slower at short times compared to the AB pair. 
This can be understood in terms of the effective poten- 
tial, discussed later. 



C. Relative diffusion: Mean square relative 
displacement (MSRD) 



fn this subsection, we investigate the time dependence 
of the mean square relative displacement (| — 
Tij(0) \ 2 )r a , the simplest physical quantity associated 
with the pair motion, where the index i and j denote 
A and/or B particles and the subscript r D indicate the 
ensemble averaging is restricted to the pairs whose ini- 
tial separation corresponds to r Q ||. First, we consider 
the case where the initial separations for the three pairs 
corresponds to the first peak of the respective partial ra- 
dial distribution functions (see figure f). In other words, 
we consider first those pairs that were initially nearest 
neighbor pairs. 

Figure 13 shows the result for the time dependence of 
the relative mean square displacement (MSRD) of the 
three pairs. At long times the MSRD varies linearly with 
time. However, the evolution of MSRD with time differs 
for different pairs. As expected, the smaller size of the B 
particles and the weak BB interaction leads to a faster 
approach of the diffusive limit of BB pair separation. The 
time scale needed to reach the diffusive limit is shorter 
for the AB pair than that for the AA pair. 

From the slope of the curves in the linear region one 
can obtain the values of the relative diffusion constants 
Dji of the different pairs. The values thus obtained 
are the following: D AA ~ 0.0032, D AB ~ 0.0048 and 
B ~ 0.0064. One should note that even though the 
difference in size of the A and B particles is small, D BB 
is almost twice of D^ A . At sufficiently long time, one 
would certainly expect the diffusion constant for the rel- 
ative motion of a pair should be the sum of the individual 
diffusion constants of the two atoms obtained from the 
slope of the corresponding mean square displacements at 
long time. Indeed, we find there is a good agreement. 



An investigation of the behavior of MSRD is also per- 
formed for atomic pairs which were initially next nearest 
neighbor. When compared to the nearest neighbor pairs 
(figure 13) we find that the slope of the corresponding 
straight lines are almost identical, although in the case 
of AA and AB pairs the diffusive limits are reached at 
shorter times. This has been shown in figure 14. One 
should remember that the AA and AB correlations are 
highest at the first coordination shell whereas the high- 
est BB correlations occur at the second coordination shell 
(see figure 1). Thus, at short time the increase in slope 
for the AA and AB pairs can be explained in terms of the 
decrease in correlations at the second coordination shell. 



D. The non-Gaussian parameter for the relative 
motion 



In a highly supercooled liquid, the single particle dis- 
placement distribution function G s (r, t) (known as the 
self-part of the van Hove correlation function) has an 
extended tail and is, in general, non-Gaussian. The 
deviation from the Gaussian behavior is often thought 
to reflect the presence of the transient inhomogencitics 
and can be characterized by the non-Gaussian parame- 
ter a 2 (i)H§ 



a 2 {t) 



3(Ar 4 (Q) 
5(Ar 2 (t)) 2 



(2) 



where (Ar 2 (t)) and (Ar (i)) are the second and fourth 
moments of G s (r, t) , respectively. At intermediate time 
scale, Q!2(t) increases with time and the maximum of 
tt2(t) occurs around the end of the /3 relaxation region. 
Only on the time scale of diffusion or the a relaxation, 
ot2{t) starts to decrease and finally at very long time limit, 
it reaches to zero, a^if) calculated for the A and B par- 
ticles are shown in figure 15. The maximum in o^W is 
seen to shift slightly towards left and also the height of 
the maximum is higher for the B particles. This is a 
clear evidence that the B particles probe much more het- 
erogeneous environment than does the A particles. This 
difference can be explained in terms of the smaller con- 
centration of B particles, different sizes of the A and B 
particles and also from the fact that the interaction be- 
tween the B particles is much weaker than that between 
the A particles (2lJ |22). 

Motivated by these findings for the single particle dis- 
placement distribution function, we introduce a new non- 
Gaussianity parameter for the pair dynamics, denoted by 
a P (t). a P (t) can be a measure of the deviation from 
the Gaussian behavior of the pair distribution function 
<72(r , r;i). It can be defined as, 



/(0) l 4 k 



5<|r«(t)-r„(0) \ 2 )i 



T, (i, j — Aand/or B) 
(3) 

where (| (t) - r^- (0) \ 2 ) To and (| r# (i) - (0) | 4 ) ro are 
the mean square relative displacement and mean quartic 
relative displacement of the ij pair. One should note that 
ce P (t) is identical to zero for a Gaussian pair distribution 
function. 

p 

In figure 16 we show the behavior of the a 2 ' 3 as a func- 
tion of time for the three different pairs. We again con- 
sider only those pairs that were initially nearest neighbor. 
The behavior observed for the three pairs is markedly dif- 
ferent. The dynamics explored by the BB pair is seen to 
be less heterogeneous than the AA and AB pairs. Be- 
cause of the smaller size of the B particles, the B par- 
ticles reach the average distribution faster, although it 
explores larger heterogeneity. The AA pair reaches the 
diffusive limit at longer time scale than that for the AB 
pair, the AB pair explore more heterogeneous dynamics 
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as is clearly evident from the difference in the maximum 
value of af (t). 



IV. THEORETICAL ANALYSIS 

For the motion of an atomic pair in a pure fluid, 
HaanQ introduced a simple mean- field level equation of 
motion for the time-dependent pair distribution function 
<72- This equation was shown to give a quantitatively cor- 
rect description both at short and long times[p^|. This 
treatment is mean-field in the sense that the two atoms 
were assumed to diffuse in an effective-force field of the 
surrounding particles given by the gradient of the poten- 
tial of mean force. The equation for gi was represented 
by a Smoluchowski equation and the correct short time 
description of 52 was obtained only by introducing a non- 
linear time that is related to the mean-squared distance 
(MSD) moved by a single atom. In other words, an ad 
hoc introduction of a time-dependent diffusion constant 
D(t) in the equation of motion gives the correct descrip- 
tion at short times. 

In the view of its success for one-component liquid, we 
have performed a similar mean-field model calculations 
for the binary mixture considered here. The generaliza- 
tion to binary mixture gives the following Smoluchowski 
equation for the different pairs 



dg2(r ,r;t) 
dm 



(r D , r; (r Q , r; t)V% (r 



( 4 ) 

where the index i and j denote A and/or B particles. (3 
is the inverse of the Boltzmann's constant (ks) times the 
absolute temperature (T). Wij(r) is the potential of mean 
force ('effective potential') between the i and j particles 
given by 

Wij(r) = -k B T\ng l3 {r) (5) 

where <?<y(r) is the partial radial distribution function. In 
Eq. 4, the "time" Ty is defined by 

' K«)-^(0)| 2 )r o 



G 



flr^-r^O) \*) + (\v j {t)-v j {Q) | 2 ) 



(G) 



where (| Yij{t) — -(0) \ 2 ) ro is the mean square relative 
displacement (MSRD) of the 'if pair. Note that an ap- 
proximation is made in the above equation by neglecting 
the cross-correlation between the two particles ('i' and 
'j') and the MSRD is replaced by the sum of individual 
particle's mean square displacement (MSD). 

Now the integration of the (r Q , r; £) over the solid 
angles Cl a and Cl corresponding to the initial and fi- 
nal positions, respectively, gives the radial part of the 
full distribution function (the zeroth-angular moment of 
g% {v ,r;t)) 

92, r ad( r o,r;t) = ^ j dn o dClgl 3 {r ,r;t) (7) 



Note that the normalization of this function is 

drr 2 gl rad {r ,r;t) = 1 (8) 





The equation of motion for 92rad( r °> r it) (derived from 
Eq. 4) is solved numerically (by Crank-Nicholson 
method) for the different pairs and the results obtained 
from this model calculations are compared with the sim- 
ulation results. The partial radial distribution functions 
gij(r) and the mean-square displacements of the A and 
B particles (required as input) are obtained from the 
present simulation. 

Figures 17 and 18 compare model calculations with 
the simulated distribution functions for the AA and AB 
nearest neighbor pairs. The time evolution of the dis- 
tributions are described well by the simple mean-field 
model. The underlying effective-potential energy sur- 
faces are plotted in figure 19. Thus relative diffusion 
in these cases can be considered as overdamped motion 
in an effective potential, which takes place mainly via 
hopping (as shown in figures 3 and 5), which governs the 
time evolution of the distributions for the AA and AB 
pair. 

Unfortunately, the good agreement observed above be- 
tween simulation and theory for the AA and the AB 
pairs, does not extend to the BB pair. This is shown 
in figure 20. As the number of B particles present in 
the system is much less (20 %) and the interparticle in- 
teraction is weak, the effective potential for a B atom 
interacting with a nearest neighbor atom is unfavorable 
(see the figure 19). Consequently, the nearest neighbor 
BB pairs execute highly anharmonic motions. Thus, the 
fluctuations about the mean-force field experienced by 
the BB pair are large and important. These fluctuations 
are neglected here, as in other mean-field level descrip- 
tion. 

The extension of the calculations to the case of next 
nearest neighbor pairs also been carried out and com- 
pared with the simulated distributions. It should be 
noted that compared to the nearest neighbor pairs, the 
AA and AB pairs are now executes motions in a rela- 
tively weak, shallow potential, whereas the motions of 
the BB pairs takes place in a relatively strong, bound 
potential well (see the figure 19). Thus, for the BB pairs, 
one expects a better agreement with the simulated dis- 
tributions compared to the earlier case (nearest neighbor 
BB pairs). Indeed, the agreement is better for the BB 
pairs, as shown in figure 21 (compare with figure 20). We 
have found that the MF model provides a good descrip- 
tion of the dynamics of the AA and AB pairs, although 
the agreement is not as satisfactory as for the nearest 
neighbor pairs. 

Thus, it is evident that the MF description for the 
time dependent pair distribution functions is reasonably 
good for the AA and AB pairs. Simulation results have 
shown that the relative diffusion of an AB pair is higher 
than that for an AA pair. We noted that this due to 
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more frequent hopping of B particles than the A par- 
ticles. Our main objective now is to see whether the 
frequent jump motions of the B particles, as predicted 
by the simulations, can be explained in terms of the MF 
model described above. 

We have performed an approximate calculation to get 
an estimate of the transition rate between the first two 
adjacent minima in the effective potential energy surface 
of the A A and AB pairs (see the figure 21). In other 
words, the rate of crossing from the deep minima located 
at the nearest neighbor pairs, to the adjacent minima 
(corresponds to the next nearest neighbors). As the mo- 
tion of a pair in the effective potential was treated by 
a Smoluchowski equation, we use the corresponding rate 
expression in the overdamped limit to calculate the es- 
cape rate. Thus, we have an expression for the escape 
rate given byp6| 



ks 



^rmn UJ max 



2ir( 



exp 



AW 
"kWf 



(9) 



where AW — W(r max ) — W(r min ) is the Arrhenius ac- 
tivation energy and Lu m i n , Lu max are the frequency at the 
minima (r TO , n ) and maxima (r max ) in the effective po- 
tential W(r), respectively. The diffusion coefficient D is 
related to the friction ( by D = ksT/(. 

Thus, to calculate the transition rate we need to know 
the values of the frequencies uj m i ni uj max , and the bar- 
rier height AW, which are different for the A A and AB 
pairs. For the AA pairs, these parameters are found to 
be Lo* min ~ 16.5, uj* max ~ 6.5 and AW A a ^ 2.25k B T, 
whereas for the AB pairs they are io* min ~ 17.8, u* max ~ 
7.4 and AW A b ^ 2.45fc s T. The relative diffusion of the 
two pairs are D AA ~ 0.0032 and D AB ~ 0.0048. Using 
all these parameters, the escape rate calculated for the 
AA and AB pairs are (in reduced units) k AA ~ 5.9 x 10~ 3 
and k AB ~ 8.8 x I0~ 3 , respectively (in terms of time r, 
which is equal to 2.2 ps for argon units). 

Even though the barrier height AWab > AWaa, the 
transition rate for the AB pair is larger than that for 
the AA pair. Thus, the jump motions are much more 
frequent for the AB pair, due to the large diffusion of the 
B particles in the potential energy surface (which mainly 
occurs via hopping mode). 



liquids in a more collective way, by following the relative 
motion of the atoms rather than absolute motion of the 
atoms. We find that the three pairs (AA, BB and AB) 
behave differently. The analysis of the trajectory shows a 
clear evidence of the jump motions for all the three pairs. 

The relative diffusion constant of the BB pair (D BB ) 
is almost twice the value for the A A pair (D^ A ). This 
clearly suggests the importance of the jump dynamics 
for the B particles and indeed, we find that the motion 
of the B particles is mostly discontinuous in nature, while 
the A particles show occasional hopping. The dynamic 
inhomogeneity present in a supercooled liquid is generally 
characterized by the well-known non-Gaussian parameter 
Q!2(t), which describe the deviations from the Gaussian 
behavior in the motion of a single atom. In this paper, 
we have generalized this concept and introduce a non- 
Gaussian parameter for the pair dynamics (a P (i)), to 
measure the deviations from the Gaussian behavior in 
the relative motion of the atoms. At intermediate times, 
all the three pair distribution functions for the three pairs 
show significant deviations from the Gaussian behavior, 
with different degree. 

It is found that for the nearest neighbor AA and AB 
pairs, which are confined to a strong effective potential 
and merely makes anharmonic motions in it, the dynam- 
ics can be treated at the mean-field level. However, as 
the motion of a nearest neighbor BB pair is highly an- 
harmonic, one must include the effects of the fluctuations 
about the mean-force field, in order to get a correct de- 
scription of the dynamics. 

While the mean-field treatment provides reasonably 
accurate description of pair dynamics (at least for AA 
and AB pairs), it must be supplemented with the time 
dependent diffusion coefficient (D(t)). This is a limita- 
tion of the mean-field approach because at present we do 
not have any theoretical means to calculate D(t) from 
first principles. The mode coupling theory (MCT) does 
not work because it neglects hopping which is the dom- 
inant mode of mass transport in deeply supercooled liq- 
uids, even when the system is quite far from the glass 
transition. As we discussed recently, the hopping can be 
coupled to anisotropy in the local stress tensor jp6[. The 
calculation of the latter is also non-trivial. Work in this 
direction is under progress. 



CONCLUSIONS 



Let us first summarize the main results of this study. 
We have presented the molecular dynamics simulation re- 
sults for the time dependent pair distribution functions in 
a strongly non-ideal glass forming binary Lennard-Jones 
mixture. In addition, a mean-field description of the pair 
dynamics is considered and the comparison is made with 
the simulated distributions. The main goal of this inves- 
tigation was to explore the dynamics of the supercooled 



Acknowledgments 

This work was supported in part by the Council of Sci- 
entific and Industrial Research (CSIR), India and the De- 
partment of Science and Technology (DST), India. One 
of the authors (R.K.M) thanks the University Grants 
Commission (UGC) for providing the Research Scholar- 
ship. 



7 



[1] I. Oppcnheim and M. Bloom, Can. J. Phys. 39, 845 
(1961). 

[2] J. A. Bucaro and T. A. Litovitz, J. Chem. Phys. 54, 3846 
(1971); A. J. C. Ladd, T. A. Litovitz, and C. J. Montrose, 
ibid. 71, 4242 (1979); A. J. C. Ladd, T. A. Litovitz, J. H. 
R. Clarke, and L. V. Woodcock, ibid. 72, 1759 (1980). 

[3] M. S. Miller, D. A. McQuarrie, G. Birnbaum, and J. D. 
Poll, J. Chem. Phys. 57, 618 (1972). 

[4] S. A. Adelman, Adv. Chem. Phys. 53, 61 (1983). 

[5] J. T. Hynes, R. Kapral, and G. M. Torrie, J. Chem. Phys. 

72, 177 (1980). 

[6] S. W. Haan, Phys. Rev. A 20, 2516 (1979). 
[7] C. D. Boley and G. Tenti, Phys. Rev. A 21, 1652 (1980). 
[8] U. Balucani and R. Vallauri, Physica A 102, 70 (1980); 
Phy. Lett. A 76, 223 (1980); Chem. Phys. Lett. 74, 75 

(1980) . 

[9] U. Balucani, R. Vallauri, C. S. Murthy, T. Gaskell, and 
M. S. Woolfson, J. Phys. C: Solid State Phys. 16, 5605 
(1983); U. Balucani, R. Vallauri, T. Gaskell, and M. Gori, 
ibid. 18, 3133 (1985). 
[10] H. Posch, F. Vesely, and W. Steele, Mol. Phys. 44, 241 

(1981) . 

[11] G. T. Evans and D. Kivelson, J. Chem. Phys. 85, 7301 
(1986). 

[12] Ten-Ming Wu and S. L. Chang, Phys. Rev. E 59, 2993 

(1999) . 

[13] G. Srinivas, A. Mukherjee, and B. Bagchi, J. Chem. Phys. 

114, 6220 (2001). 
[14] R. K. Murarka and B. Bagchi, J. Chem. Phys. 117, 1155 

(2002). 

[15] J. G. Saven and J. L. Skinner, J. Chem. Phys. 99, 4391 
(1993); M. D. Stephens, J. G. Saven, and J. L. Skinner, 
J. Chem. Phys. 106, 2129 (1997). 

[16] A. D. Santis, A. Ercoli, and D. Rocca, Phys. Rev. E 57, 
R4871 (1998); Phys. Rev. Lett. 82, 3452 (1999). 

[17] M. D. Ediger, C. A. Angell, and S. R. Nagel, J. Phys. 
Chem. 100, 13200 (1996); C. A. Angell, Science 267, 
1924 (1995); C. A. Angell, K. L. Ngai, G. B. McKenna, P. 
F. McMillan, and S. W. Martin, J. Appl. Phys. 88, 3113 

(2000) ; P. G. Debenedetti and F. H. Stillinger, Nature 
410, 259 (2001). 

[18] A. Mukherjee, S. Bhattacharyya, and B. Bagchi, J. 
Chem. Phys. 116, 4577 (2002). 

[19] M. D. Ediger, Annu. Rev. Phys. Chem. 51, 99 (2000). 

[20] E. R. Weeks, J. C. Crocker, A. C. Levitt, A. Schofield, 
and D. A. Weitz, Science 287, 627 (2000); W. K. Kegel 
and A. van Blaaderen, Science 287, 290 (2000); L. A. De- 
schenes and D. A. Vanden Bout, Science 292, 255 (2001). 

[21] W. Kob and H. C. Andersen, Phys. Rev. E 51, 4626 
(1995); W. Kob and H. C. Andersen, Phys. Rev. Lett. 

73, 1376 (1994). 

[22] W. Kob, C. Donati, S. J. Plimpton, P. H. Poole, and S. 
C. Glotzer, Phys. Rev. Lett. 79, 2827 (1997); C. Donati, 
S. C. Glotzer, P. H. Poole, W. Kob, and S. J. Plimpton, 
Phys. Rev. E 60, 3107 (1999); K. Vollmayr-Lee, W. Kob, 
K. Binder, and A. Zippelius, J. Chem. Phys. 116, 5158 
(2002). 

[23] S. Sastry, P. G. Debenedetti, and F. H. Stillinger, Nature 

393, 554 (1998). 
[24] B. Doliwa and A. Heuer, Phys. Rev. Lett. 80, 4915 

(1998); J. Qian, R. Hentschke, and A. Heuer, J. Chem. 



Phys. 110, 4514 (1999). 
[25] C. Donati, J. F. Douglas, W. Kob, S. J. Plimpton, P. 
H. Poole, and S. C. Glotzer, Phys. Rev. Lett. 80, 2338 

(1998) . 

[26] S. Bhattacharyya and B. Bagchi, Phys. Rev. Lett. 89, 
025504 (2002). 

[27] S. Bhattacharyya, A. Mukherjee, and B. Bagchi, J. 

Chem. Phys. 117, 2741 (2002). 
[28] C. Oligschleger and H. R. Schober, Phys. Rev. B 59, 811 

(1999) . 

[29] H. Miyagawa, Y. Hiwatari, B. Bernu, and J. P. Hansen, 

J. Chem. Phys. 88, 3879 (1988). 
[30] D. Caprion, J. Matsui, and H. R. Schober, Phys. Rev. 

Lett. 85, 4293 (2000). 
[31] M. P. Allen and D. J. Tildesley, Computer Simulation of 

Liquids (Oxford University Press, Oxford, 1987). 
[32] G. J. Martyna, D. J. Tobias, and M. L. Klein, J. Chem. 

Phys. 101, 4177 (1994); H. C. Andersen, J. Chem. Phys. 

72, 2384 (1980); S. Nose, Mol. Phys. 52, 255 (1984); W. 

G. Hoover, Phys. Rev. A 31, 1695 (1985). 
[33] G. J. Martyna, M. E. Tuckerman, and M. L. Klein, J. 

Chem. Phys. 97, 2635 (1992). 
[34] M. E. Tuckerman, G. J. Martyna, and B. J. Berne, J. 

Chem. Phys. 97, 1990 (1992). 
[35] G. J. Martyna, M. E. Tuckerman, D. J. Tobias, and M. 

L. Klein, Mol. Phys. 87, 1117 (1996). 
[36] R. Zwanzig, N on- equilibrium Statistical Mechanics (Ox- 
ford University Press, New York, 2001). 



8 



Figure Captions 

Figure 1. The radial distribution function g(r) for the 
AA, AB, and BB correlation is plotted against distance. 
The solid line is gAA(r), the dashed line is gAB {r), and 
the dot-dashed line is .gssM- For details, see the text. 

Figure 2. The radial part of the time-dependent pair 
distribution function g2,rad{r 0l 'i';t) for the AA pair as 
a function of separation r at four different times: (a) 
t = 20r, (b) t = 50t, (c) t = lOOr, and (d) t = 300r. The 
initial separation r a corresponds to the first maximum of 
gAA{r)- Note that the time unit r = 2.2ps if argon units 
are assumed. 

Figure 3. Projections into x-y plane of the trajectory 
of a typical nearest neighbor AA pair over a time interval 
t = 500t. Note that the time unit r = 2.2ps for argon 
units. 

Figure 4. The radial part of the time-dependent pair 
distribution function g2,rad(r ,i r ;t) for the AB pair as 
a function of separation r at four different times: (a) 
t = 20r, (b) t = 50r, (c) t = lOOr, and (d) t = 300r. The 
initial separation r corresponds to the first maximum of 
gAB{r)- The time unit r = 2.2ps for the argon units. For 
further details, see the text. 

Figure 5. (a) Projections into x-y plane of the tra- 
jectory of a typical nearest neighbor AB pair over a time 
interval t = 500t. (b) Trajectory of the A particle of the 
same AB pair as in (a), within the same time window, 
(c) Trajectory of the B particle of the same AB pair. 

Figure 6. The net displacement of an AB pair into 
x-y plane (AL xy ) as shown in figure 5(a), in the same 
time interval. Note that the displacement is quite large. 

Figure 7. The radial part of the time-dependent pair 
distribution function g2.rad{i'o 1 r;t) for the BB pair as 
a function of separation r at four different times: (a) 
t = 20t, (b) t = 50r, (c) t = lOOr, and (d) t = 300r. 
The initial separation r D corresponds to the first peak of 
9bb{t)- 

Figure 8. Projections into x-y plane of the trajectory 
of a typical nearest neighbor BB pair over a time interval 
t = 500t. 

Figure 9. The radial part of the pair distribution 
function g2,rad(r ,r;t) for the AA pair at four different 
times: (a) t = 4r, (b) t = 20r, (c) t = lOOr, and (d) 
t = 300t. Here the initial separation r Q is chosen at the 
second peak of gAA^r)- 

Figure 10. The radial part of the pair distribution 
function g2,rad(? ,r;t) for the AB pair at four different 
times: (a) t = At, (b) t = 20t, (c) t = lOOr, and (d) 
t = 300t. Here the initial separation r a corresponds to 
the second peak of gAB{r)- 

Figure 11. The radial part of the pair distribution 
function g2,rad(r , r; t) for the BB pair at four different 
times: (a) t = At, (b) t = 20t, (c) t = lOOr, and (d) 
t = 300r. Here the initial separation r D is chosen at the 
second peak of gBB{r)- 

Figure 12. (a) The angular part of the time- 
dependent pair distribution function g2,ang 

(r ,6;t) for 

the AA pair at four different times, (b) g2,ang{r ,6;t) 



for the AB pair, (c) g2,ang{r , 0; t) for the BB pair. In all 
the three cases, we consider only those pairs which were 
initially separated at the nearest neighbor distance. For 
further details, see the text. 

Figure 13. Time dependence of the mean square rel- 
ative displacement (MSRD) for the AA, AB and BB pair 
(in units of (Jaa)- The initial separation r of the three 
pairs corresponds to the first peak of the respective par- 
tial radial distribution functions. The solid line repre- 
sents the result for the AA pair, the dashed line AB pair, 
and the dotted line for the BB pair. For the detailed 
discussion, see the text. 

Figure 14. (a) Comparison of the mean square rela- 
tive displacement (MSRD) for the AA pair with different 
initial separations. The solid line represents the nearest 
neighbor AA pair and the dashed line represents the next 
nearest neighbor AA pair, (b) Same as in (a), but for the 
AB pair, (c) For the BB pair. For details, see the text. 

Figure 15. The behavior of the non-Gaussian param- 
eter ct2 (t) as a function of time for the A and B particles. 
The solid line is for the A particles and the dashed line 
for the B particles. 

Figure 16. The behavior of the non-Gaussian param- 
eter a f (t) as a function of time for the AA, AB, and BB 
pairs, initially separated at the nearest neighbor distance. 
The solid line represents the result for the AA pair, the 
dashed line for the AB pair, and the dot-dashed line for 
the BB pair. 

Figure 17. The simulated distribution g2,rad(r ,r;t) 
for the AA pair is compared with the mean-field model 
calculations at three different times: (a) t = I0r, (b) 
t = 50t, and (c) t = lOOr. The initial separation r a cor- 
responds to the first maximum of gAA {f) ■ The histogram 
represents the simulation results and the dashed line rep- 
resents the results of the model calculations. Note that 
the time unit r = 2.2ps if argon units are assumed. 

Figure 18. The simulated distribution g2,rad{r ,r;t) 
for the AB pair is compared with the mean-field model 
calculations at three different times: (a) t = lOr, (b) 
t = 50t, and (c) t = lOOr. The initial separation r Q 
corresponds to the first maximum of gAB (r). The his- 
togram represents the simulation results and the dashed 
line represents the results of the model calculations. 

Figure 19. The potential of mean force W(r) for the 
AA, AB and BB pairs in the Kob- Andersen model at the 
reduced pressure P* = 20 and the reduced temperature 
T* = 1.0. The solid line represents for the AA pair, the 
dashed line for the AB pair, and the dot-dashed line for 
the BB pair. 

Figure 20. The simulated distribution g2,rad{r ,r;t) 
for the BB pair is compared with the mean-field model 
calculations at three different times: (a) t = lOr, (b) 
t = 50t, and (c) t = lOOr. The initial separation r Q 
corresponds to the first peak of 5b_b(^)- The histogram 
represents the simulation results and the dashed line rep- 
resents the results of the model calculations. 

Figure 21. The simulated distribution g2,rad(r ,'r;t) 
for the BB pair is compared with the mean-field model 
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calculations at two different times: (a) t — lOrand (b) sents the simulation results and the dashed line repre- 
t = IOOt. Here the initial separation r a corresponds to sents the results of the model calculations. For further 
the second peak of SfBs(r). The histogram again repre- details, see the text. 
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